######################################
######################################
######################################
###Robustness Models -- All Disease###
######################################
######################################
######################################

###########################
####Affected Individuals###
###########################
dis.dat$Affected1000 <- dis.dat$Affected/1000
dis.dat$lagAffected1000 <- dis.dat$lagAffected/1000
####Cattle
iv2.cat.aff <-  ivreg(Affected1000~cattle_meat_tonnes1000000+t+
                    lagAffected1000+logpopulation+logGDPpc+factor(ccode)|
                    comb.np.use_tonnes1000000+t+lagAffected1000+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.aff, diagnostics=T)
#Cluster results
iv2.cat.aff.sum <- coeftest(iv2.cat.aff, vcov = vcovCL(iv2.cat.aff, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("Affected1000","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagAffected1000","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagAffected1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagAffected1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.aff.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.aff <-  ivreg(Affected1000~chicken_meat_tonnes1000000+t+
                      lagAffected1000+logpopulation+logGDPpc+factor(ccode)|
                      comb.np.use_tonnes1000000+t+lagAffected1000+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.aff, diagnostics=T)
#Cluster results
iv2.chick.aff.sum <- coeftest(iv2.chick.aff, vcov = vcovCL(iv2.chick.aff, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("Affected1000","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagAffected1000","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagAffected1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagAffected1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.aff.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.aff <-  ivreg(Affected1000~pigs_meat_tonnes1000000+t+lagAffected1000+logpopulation+logGDPpc+factor(ccode)|
                     comb.np.use_tonnes1000000+t+lagAffected1000+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.aff, diagnostics=T)
#Cluster results
iv2.pigs.aff.sum <- coeftest(iv2.pigs.aff, vcov = vcovCL(iv2.pigs.aff, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("Affected1000","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagAffected1000","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagAffected1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagAffected1000+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.aff.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.aff.sum, iv2.chick.aff.sum, iv2.pigs.aff.sum)
#Obs and diagnostics
stargazer(iv2.cat.aff, iv2.chick.aff, iv2.pigs.aff)
#Corrected weak instruments
iv.main.aff <- data.frame(cbind(iv2.cat.aff.sum.stock,iv2.chick.aff.sum.stock,iv2.pigs.aff.sum.stock))
colnames(iv.main.aff) <- c("Cattle (Affected)", "Chicken (Affected)", "Pigs (Affected)")
iv.main.aff


##############################
###Cause (Poultry and Pigs)###
##############################
####Chicken
iv2.chick.cause <-  ivreg(cause_poultry~chicken_meat_tonnes1000000+t+
                            lagcause_poultry+logpopulation+logGDPpc+factor(ccode)|
                            comb.np.use_tonnes1000000+t+lagcause_poultry+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.cause, diagnostics=T)
#Cluster results
iv2.chick.cause.sum <- coeftest(iv2.chick.cause, vcov = vcovCL(iv2.chick.cause, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("cause_poultry","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagcause_poultry","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagcause_poultry+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagcause_poultry+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.cause.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.cause <-  ivreg(cause_pig~pigs_meat_tonnes1000000+t+lagcause_pig+logpopulation+logGDPpc+factor(ccode)|
                           comb.np.use_tonnes1000000+t+lagcause_pig+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.cause, diagnostics=T)
#Cluster results
iv2.pigs.cause.sum <- coeftest(iv2.pigs.cause, vcov = vcovCL(iv2.pigs.cause, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("cause_pig","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagcause_pig","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagcause_pig+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagcause_pig+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.cause.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.chick.cause.sum, iv2.pigs.cause.sum)
#Obs and diagnostics
stargazer(iv2.chick.cause, iv2.pigs.cause)
#Corrected weak instruments
iv.main.cause <- data.frame(cbind(iv2.chick.cause.sum.stock,iv2.pigs.cause.sum.stock))
colnames(iv.main.cause) <- c("Chicken (Cause)", "Pigs (Cause)")
iv.main.cause


#######################################
#######################################
###Alternative Explanatory Variables###
#######################################
#######################################

###########################
###Meat Production Index###
###########################
#Change Bangladesh and Brunei (Muslim states) to 0
dis.dat$pigs_meat_PI <- ifelse((dis.dat$ccode==771|dis.dat$ccode==835), 0, dis.dat$pigs_meat_PI)
####Cattle
iv2.cat.PI <-  ivreg(outbreak_event~cattle_meat_PI+t+
                    lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                    comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.PI, diagnostics=T)
#Cluster results
iv2.cat.PI.sum <- coeftest(iv2.cat.PI, vcov = vcovCL(iv2.cat.PI, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_PI","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(cattle_meat_PI ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_PI ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.PI.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.PI <-  ivreg(outbreak_event~chicken_meat_PI+t+
                      lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                      comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.PI, diagnostics=T)
#Cluster results
iv2.chick.PI.sum <- coeftest(iv2.chick.PI, vcov = vcovCL(iv2.chick.PI, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_PI","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(chicken_meat_PI ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_PI ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.PI.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.PI <-  ivreg(outbreak_event~pigs_meat_PI+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                     comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.PI, diagnostics=T)
#Cluster results
iv2.pigs.PI.sum <- coeftest(iv2.pigs.PI, vcov = vcovCL(iv2.pigs.PI, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","pigs_meat_PI","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(pigs_meat_PI ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_PI ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.PI.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.PI.sum, iv2.chick.PI.sum, iv2.pigs.PI.sum)
#Obs and diagnostics
stargazer(iv2.cat.PI, iv2.chick.PI, iv2.pigs.PI)
#Corrected weak instruments
iv.main.PI <- data.frame(cbind(iv2.cat.PI.sum.stock,iv2.chick.PI.sum.stock,iv2.pigs.PI.sum.stock))
colnames(iv.main.PI) <- c("Cattle (PI)", "Chicken (PI)", "Pigs (PI)")
iv.main.PI

#########################################
#########################################
###Disaggregate Instrumental Variables###
#########################################
#########################################


##############
###Nitrates###
##############
####Cattle
iv2.cat.nit <-  ivreg(outbreak_event~cattle_meat_tonnes1000000+t+
                    lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                    nitrogen.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.nit, diagnostics=T)
#Cluster results
iv2.cat.nit.sum <- coeftest(iv2.cat.nit, vcov = vcovCL(iv2.cat.nit, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","t","ccode","nitrogen.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ nitrogen.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.nit.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.nit <-  ivreg(outbreak_event~chicken_meat_tonnes1000000+t+
                      lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                      nitrogen.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.nit, diagnostics=T)
#Cluster results
iv2.chick.nit.sum <- coeftest(iv2.chick.nit, vcov = vcovCL(iv2.chick.nit, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_tonnes1000000","t","ccode","nitrogen.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ nitrogen.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.nit.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.nit <-  ivreg(outbreak_event~pigs_meat_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                     nitrogen.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.nit, diagnostics=T)
#Cluster results
iv2.pigs.nit.sum <- coeftest(iv2.pigs.nit, vcov = vcovCL(iv2.pigs.nit, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","pigs_meat_tonnes1000000","t","ccode","nitrogen.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ nitrogen.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.nit.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.nit.sum, iv2.chick.nit.sum, iv2.pigs.nit.sum)
#Obs and diagnostics
stargazer(iv2.cat.nit, iv2.chick.nit, iv2.pigs.nit)
#Corrected weak instruments
iv.main.nit <- data.frame(cbind(iv2.cat.nit.sum.stock,iv2.chick.nit.sum.stock,iv2.pigs.nit.sum.stock))
colnames(iv.main.nit) <- c("Cattle (Nitrate)", "Chicken (Nitrate)", "Pigs (Nitrate)")
iv.main.nit

################
###Phosphates###
################

####Cattle
iv2.cat.phos <-  ivreg(outbreak_event~cattle_meat_tonnes1000000+t+
                    lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                    phosphate.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.phos, diagnostics=T)
#Cluster results
iv2.cat.phos.sum <- coeftest(iv2.cat.phos, vcov = vcovCL(iv2.cat.phos, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","t","ccode","phosphate.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ phosphate.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.phos.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.phos <-  ivreg(outbreak_event~chicken_meat_tonnes1000000+t+
                      lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                      phosphate.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.phos, diagnostics=T)
#Cluster results
iv2.chick.phos.sum <- coeftest(iv2.chick.phos, vcov = vcovCL(iv2.chick.phos, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_tonnes1000000","t","ccode","phosphate.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ phosphate.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.phos.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.phos <-  ivreg(outbreak_event~pigs_meat_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                     phosphate.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.phos, diagnostics=T)
#Cluster results
iv2.pigs.phos.sum <- coeftest(iv2.pigs.phos, vcov = vcovCL(iv2.pigs.phos, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","pigs_meat_tonnes1000000","t","ccode","phosphate.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ phosphate.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.phos.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.phos.sum, iv2.chick.phos.sum, iv2.pigs.phos.sum)
#Obs and diagnostics
stargazer(iv2.cat.phos, iv2.chick.phos, iv2.pigs.phos)
#Corrected weak instruments
iv.main.phos <- data.frame(cbind(iv2.cat.phos.sum.stock,iv2.chick.phos.sum.stock,iv2.pigs.phos.sum.stock))
colnames(iv.main.phos) <- c("Cattle (Phosphates)", "Chicken (Phosphates)", "Pigs (Phosphates)")
iv.main.phos


##############################
##############################
###Controls and Confounders###
##############################
##############################


#########################
####Imports and Export###
#########################

##Create log values
dis.dat$logcattle.meat.imp_1kusd <- log(dis.dat$cattle.meat.imp_1kusd+1)
dis.dat$logchicken.meat.imp_1kusd <- log(dis.dat$chicken.meat.imp_1kusd+1)
dis.dat$logpigs.meat.imp_1kusd <- log(dis.dat$pigs.meat.imp_1kusd+1)
dis.dat$logcattle.meat.exp_1kusd <- log(dis.dat$cattle.meat.exp_1kusd+1)
dis.dat$logchicken.meat.exp_1kusd <- log(dis.dat$chicken.meat.exp_1kusd+1)
dis.dat$logpigs.meat.exp_1kusd <- log(dis.dat$pigs.meat.exp_1kusd+1)


####Cattle
iv2.cat.impexp <-  ivreg(outbreak_event~cattle_meat_tonnes1000000+t+
                    lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+
                    logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                    logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd|
                    comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+
                    logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                    logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data=dis.dat)
summary(iv2.cat.impexp, diagnostics=T)
#Cluster results
iv2.cat.impexp.sum <- coeftest(iv2.cat.impexp, vcov = vcovCL(iv2.cat.impexp, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc", "logcattle.meat.exp_1kusd","logchicken.meat.exp_1kusd","logpigs.meat.exp_1kusd","logcattle.meat.imp_1kusd","logchicken.meat.imp_1kusd","logpigs.meat.imp_1kusd")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.impexp.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.impexp <-  ivreg(outbreak_event~chicken_meat_tonnes1000000+t+
                      lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                      logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd|
                      comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+
                      logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                      logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data=dis.dat)
summary(iv2.chick.impexp, diagnostics=T)
#Cluster results
iv2.chick.impexp.sum <- coeftest(iv2.chick.impexp, vcov = vcovCL(iv2.chick.impexp, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc","logcattle.meat.exp_1kusd","logchicken.meat.exp_1kusd","logpigs.meat.exp_1kusd","logcattle.meat.imp_1kusd","logchicken.meat.imp_1kusd","logpigs.meat.imp_1kusd")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.impexp.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.impexp <-  ivreg(outbreak_event~pigs_meat_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                     logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd|
                     comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+
                     logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
                     logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data=dis.dat)
summary(iv2.pigs.impexp, diagnostics=T)
#Cluster results
iv2.pigs.impexp.sum <- coeftest(iv2.pigs.impexp, vcov = vcovCL(iv2.pigs.impexp, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc","logcattle.meat.exp_1kusd","logchicken.meat.exp_1kusd","logpigs.meat.exp_1kusd","logcattle.meat.imp_1kusd","logchicken.meat.imp_1kusd","logpigs.meat.imp_1kusd")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)+
          logcattle.meat.exp_1kusd+logchicken.meat.exp_1kusd+logpigs.meat.exp_1kusd+
          logcattle.meat.imp_1kusd+logchicken.meat.imp_1kusd+logpigs.meat.imp_1kusd, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.impexp.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.impexp.sum, iv2.chick.impexp.sum, iv2.pigs.impexp.sum)
#Obs and diagnostics
stargazer(iv2.cat.impexp, iv2.chick.impexp, iv2.pigs.impexp)
#Corrected weak instruments
iv.main.impexp <- data.frame(cbind(iv2.cat.impexp.sum.stock,iv2.chick.impexp.sum.stock,iv2.pigs.impexp.sum.stock))
colnames(iv.main.impexp) <- c("Cattle (Imp./Exp.)", "Chicken (Imp./Exp.)", "Pigs (Imp./Exp.)")
iv.main.impexp

############################
###Local Consumption Only###
############################
#Normalized by 1M
dis.dat$cattle.meat.exp.ton1000000 <- (dis.dat$cattle.meat.exp.ton/1000000)
dis.dat$chicken.meat.exp.ton1000000 <- (dis.dat$chicken.meat.exp.ton/1000000)
dis.dat$pigs.meat.exp.ton1000000 <- (dis.dat$pigs.meat.exp.ton/1000000)
#Create local consumption indicators
dis.dat$cattle_meat_tonnes1000000_local <- dis.dat$cattle_meat_tonnes1000000-dis.dat$cattle.meat.exp.ton1000000
dis.dat$chicken_meat_tonnes1000000_local <- dis.dat$chicken_meat_tonnes1000000-dis.dat$chicken.meat.exp.ton1000000
dis.dat$pigs_meat_tonnes1000000_local <- dis.dat$pigs_meat_tonnes1000000-dis.dat$pigs.meat.exp.ton1000000
#Correct negative error
dis.dat$cattle_meat_tonnes1000000_local <- ifelse(dis.dat$cattle_meat_tonnes1000000_local<0,0,dis.dat$cattle_meat_tonnes1000000_local)
summary(dis.dat$cattle_meat_tonnes1000000_local)
summary(dis.dat$chicken_meat_tonnes1000000_local)
summary(dis.dat$pigs_meat_tonnes1000000_local)

####Run models
####Cattle
iv2.cat.local <-  ivreg(outbreak_event~cattle_meat_tonnes1000000_local+t+
                          lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                          comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.local, diagnostics=T)
#Cluster results
iv2.cat.local.sum <- coeftest(iv2.cat.local, vcov = vcovCL(iv2.cat.local, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000_local","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000_local ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000_local ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.local.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.local <-  ivreg(outbreak_event~chicken_meat_tonnes1000000_local+t+
                            lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                            comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.local, diagnostics=T)
#Cluster results
iv2.chick.local.sum <- coeftest(iv2.chick.local, vcov = vcovCL(iv2.chick.local, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_tonnes1000000_local","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000_local ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000_local ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.local.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.local <-  ivreg(outbreak_event~pigs_meat_tonnes1000000_local+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                           comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.local, diagnostics=T)
#Cluster results
iv2.pigs.local.sum <- coeftest(iv2.pigs.local, vcov = vcovCL(iv2.pigs.local, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","pigs_meat_tonnes1000000_local","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000_local ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000_local ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.local.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.local.sum, iv2.chick.local.sum, iv2.pigs.local.sum)
#Obs and diagnostics
stargazer( iv2.cat.local, iv2.chick.local, iv2.pigs.local)
#Corrected weak instruments
iv.main.local <- data.frame(cbind(iv2.cat.local.sum.stock,iv2.chick.local.sum.stock,iv2.pigs.local.sum.stock))
colnames(iv.main.local) <- c("Cattle (Local)", "Chicken (Local)", "Pigs (Local)")
iv.main.local

#####################
####No Country FEs###
#####################
####Cattle
iv2.cat.nfes <-  ivreg(outbreak_event~cattle_meat_tonnes1000000+t+
                    lagoutbreak_event+logpopulation+logGDPpc|
                    comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc, data=dis.dat)
summary(iv2.cat.nfes, diagnostics=T)
#Cluster results
iv2.cat.nfes.sum <- coeftest(iv2.cat.nfes, vcov = vcovCL(iv2.cat.nfes, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.nfes.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.nfes <-  ivreg(outbreak_event~chicken_meat_tonnes1000000+t+
                      lagoutbreak_event+logpopulation+logGDPpc|
                      comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc, data=dis.dat)
summary(iv2.chick.nfes, diagnostics=T)
#Cluster results
iv2.chick.nfes.sum <- coeftest(iv2.chick.nfes, vcov = vcovCL(iv2.chick.nfes, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.nfes.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.nfes <-  ivreg(outbreak_event~pigs_meat_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc|
                     comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc, data=dis.dat)
summary(iv2.pigs.nfes, diagnostics=T)
#Cluster results
iv2.pigs.nfes.sum <- coeftest(iv2.pigs.nfes, vcov = vcovCL(iv2.pigs.nfes, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc, data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc, data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.nfes.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.nfes.sum, iv2.chick.nfes.sum, iv2.pigs.nfes.sum)
#Obs and diagnostics
stargazer(iv2.cat.nfes, iv2.chick.nfes, iv2.pigs.nfes)
#Corrected weak instruments
iv.main.nfes <- data.frame(cbind(iv2.cat.nfes.sum.stock,iv2.chick.nfes.sum.stock,iv2.pigs.nfes.sum.stock))
colnames(iv.main.nfes) <- c("Cattle (NFEs)", "Chicken (NFEs)", "Pigs (NFEs)")
iv.main.nfes


####################
####Year Only FEs###
####################
####Cattle
iv2.cat.yfe <-  ivreg(outbreak_event~cattle_meat_tonnes1000000+
                    lagoutbreak_event+logpopulation+logGDPpc+factor(year)|
                    comb.np.use_tonnes1000000+lagoutbreak_event+logpopulation+logGDPpc+factor(year), data=dis.dat)
summary(iv2.cat.yfe, diagnostics=T)
#Cluster results
iv2.cat.sum.yfe <- coeftest(iv2.cat.yfe, vcov = vcovCL(iv2.cat.yfe, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc","year")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+lagoutbreak_event+logpopulation+logGDPpc+factor(year), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ lagoutbreak_event+logpopulation+logGDPpc+factor(year), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.sum.yfe.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.yfe <-  ivreg(outbreak_event~chicken_meat_tonnes1000000+
                      lagoutbreak_event+logpopulation+logGDPpc+factor(year)|
                      comb.np.use_tonnes1000000+lagoutbreak_event+logpopulation+logGDPpc+factor(year), data=dis.dat)
summary(iv2.chick.yfe, diagnostics=T)
#Cluster results
iv2.chick.sum.yfe <- coeftest(iv2.chick.yfe, vcov = vcovCL(iv2.chick.yfe, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc","year")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+lagoutbreak_event+logpopulation+logGDPpc+factor(year), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ lagoutbreak_event+logpopulation+logGDPpc+factor(year), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.sum.yfe.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.yfe <-  ivreg(outbreak_event~pigs_meat_tonnes1000000+
                         lagoutbreak_event+logpopulation+logGDPpc+factor(year)|
                     comb.np.use_tonnes1000000+lagoutbreak_event+logpopulation+logGDPpc+factor(year), data=dis.dat)
summary(iv2.pigs.yfe, diagnostics=T)
#Cluster results
iv2.pigs.sum.yfe <- coeftest(iv2.pigs.yfe, vcov = vcovCL(iv2.pigs.yfe, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc","year")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+lagoutbreak_event+logpopulation+logGDPpc+factor(year), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ lagoutbreak_event+logpopulation+logGDPpc+factor(year), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.sum.yfe.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.sum.yfe, iv2.chick.sum.yfe, iv2.pigs.sum.yfe)
#Obs and diagnostics
stargazer(iv2.cat.yfe, iv2.chick.yfe, iv2.pigs.yfe)
#Corrected weak instruments
iv.main.nfes.yfe <- data.frame(cbind(iv2.cat.sum.yfe.stock,iv2.chick.sum.yfe.stock,iv2.pigs.sum.yfe.stock))
colnames(iv.main.nfes.yfe) <- c("Cattle (YFEs)", "Chicken (YFEs)", "Pigs (YFEs)")
iv.main.nfes.yfe

###########################
####Country and Year FEs###
###########################
####Cattle
iv2.cat.cyf <-  ivreg(outbreak_event~cattle_meat_tonnes1000000+factor(year)+
                    logpopulation+logGDPpc+factor(ccode)|
                    comb.np.use_tonnes1000000+factor(year)+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.cyf, diagnostics=T)
#Cluster results
iv2.cat.sum.cyf <- coeftest(iv2.cat.cyf, vcov = vcovCL(iv2.cat.cyf, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc","year")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+factor(year)+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ factor(year)+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.sum.stock.cyf <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.cyf <-  ivreg(outbreak_event~chicken_meat_tonnes1000000+factor(year)+
                      logpopulation+logGDPpc+factor(ccode)|
                      comb.np.use_tonnes1000000+factor(year)+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.cyf, diagnostics=T)
#Cluster results
iv2.chick.sum.cyf <- coeftest(iv2.chick.cyf, vcov = vcovCL(iv2.chick.cyf, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc","year")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+factor(year)+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ factor(year)+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.sum.stock.cyf <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.cyf <-  ivreg(outbreak_event~pigs_meat_tonnes1000000+factor(year)+logpopulation+logGDPpc+factor(ccode)|
                     comb.np.use_tonnes1000000+factor(year)+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.cyf, diagnostics=T)
#Cluster results
iv2.pigs.sum.cyf <- coeftest(iv2.pigs.cyf, vcov = vcovCL(iv2.pigs.cyf, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc","year")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+factor(year)+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ factor(year)+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.sum.stock.cyf <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.sum.cyf, iv2.chick.sum.cyf, iv2.pigs.sum.cyf)
#Obs and diagnostics
stargazer(iv2.cat.cyf, iv2.chick.cyf, iv2.pigs.cyf)
#Corrected weak instruments
iv.main.nfes.cyf <- data.frame(cbind(iv2.cat.sum.stock.cyf,iv2.chick.sum.stock.cyf,iv2.pigs.sum.stock.cyf))
colnames(iv.main.nfes.cyf) <- c("Cattle (CYFEs)", "Chicken (CYFEs)", "Pigs (CYFEs)")
iv.main.nfes.cyf

################
################
###GMM Models###
################
################
library(plm)

##Keep only variables of interest
subset.dat1 <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","lagoutbreak_event","logpopulation","logGDPpc","ccode", "year")]) 
#Convert data to pgmm formal
gmm.dat <- pdata.frame(subset.dat1, index=c("ccode", "year"))


#set seed
set.seed(50)

###Cattle
gmm.cattle <- pgmm(outbreak_event~cattle_meat_tonnes1000000+t+
                     lagoutbreak_event+logpopulation+logGDPpc| 
                     lag(outbreak_event,2:99),
                   data = gmm.dat, 
                   effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.cattle)

###Cattle
gmm.chicken <- pgmm(outbreak_event~chicken_meat_tonnes1000000+t+
                      lagoutbreak_event+logpopulation+logGDPpc| 
                      lag(outbreak_event,2:99),
                    data = gmm.dat, 
                    effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.chicken)

###Pigs
gmm.pigs <- pgmm(outbreak_event~pigs_meat_tonnes1000000+t+
                   lagoutbreak_event+logpopulation+logGDPpc| 
                   lag(outbreak_event,2:99),
                 data = gmm.dat, 
                 effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.pigs)

####Export to LaTex
stargazer(gmm.cattle, gmm.chicken,gmm.pigs)






